% Copyright (C) 2021-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% For a copy of the GNU General Public License,
% see <http://www.gnu.org/licenses/>.

%% Generates Figure B.2: Monte Carlo irf – full employment

% Display moments of float and peg
addpath('Auxiliary_functions/')
    
% clear;
loaded_IRFs.peg_ergodic=load('IRF_results/irfs_peg_ergodic_MC.mat');
ending_cell={'peg_ergodic','peg'};

rng(2);
loaded_PFI_peg=load('Policy_functions/pfi_peg_g_low_beta_large_grid.mat');

T_ergodic = 4e4; %periods for simulation of mean;
T_burnin = 1e4; %burn-in period

IRF_length=9;

if ~isfield(loaded_PFI_peg.par,'phi_e')
    loaded_PFI_peg.par.phi_e= 0; %exchange rate policy parameter
end

%% Set initial conditions of exogenous states, debt, and past wage equal to their respective unconditional means
temp=unique(loaded_PFI_peg.grid.y_level);
y_0 = temp(round(length(temp)/2));
temp=unique(loaded_PFI_peg.grid.r_level);
r_0 = temp(round(length(temp)/2));
temp=unique(loaded_PFI_peg.grid.g_level);
g_0 = temp(round(length(temp)/2));
d_0=loaded_PFI_peg.grid.d_grid(round(length(loaded_PFI_peg.grid.d_grid)/2));
p_0=2;
w_0=loaded_PFI_peg.grid.w_grid_vec(round(length(loaded_PFI_peg.grid.w_grid_vec)/2));

%P_trans is the transition probability matrix of tradable output. Cumul_P_trans is the cumulative probability matrix (useful for drawing tradable output realizations)
Cum_prob = cumsum(loaded_PFI_peg.MC.P_trans,2);

[Y_base_peg, var_name, pos]=simulate_model_general(y_0,r_0,g_0,d_0,w_0,p_0,T_ergodic,loaded_PFI_peg.dp,loaded_PFI_peg.h,loaded_PFI_peg.grid,loaded_PFI_peg.par,Cum_prob,1);
Y_base_peg=Y_base_peg(T_burnin+1:end,:);


G=Y_base_peg(:,pos.g_N_rule);
G_base=Y_base_peg(:,pos.g_N);
Y=Y_base_peg(:,pos.y_aggregate);
Y_N=Y_base_peg(:,pos.y_N);
RER=Y_base_peg(:,pos.rer);
W=Y_base_peg(:,pos.w);
U=Y_base_peg(:,pos.unemployment);

case_iter=1;

switch case_iter
    case 1
        [b.G_shock,bint,G_shock,rint,stats]=regress((G(2:end)),[ones(length(G(2:end)),1), (G_base(1:end-1)), Y(1:end-1), W(1:end-1)]);
    case 2
        [b.G_shock,bint,G_shock,rint,stats]=regress(log(G(2:end)),[ones(length(G(2:end)),1), G_base(1:end-1), RER(1:end-1)]);
end

G_positive=zeros(size(G_shock));
G_positive(G_shock>0)=G_shock(G_shock>0);
G_negative=zeros(size(G_shock));
G_negative(G_shock<0)=G_shock(G_shock<0);

cut.G=G(2:end);
cut.Y=Y(2:end);
cut.Y_N=Y_N(2:end);
cut.RER=RER(2:end);
cut.U=U(2:end);
b.coeff.Y=NaN(IRF_length,2);
b.se.Y=NaN(IRF_length,2);
b.coeff.G=NaN(IRF_length,2);
b.se.G=NaN(IRF_length,2);
b.coeff.RER=NaN(IRF_length,2);
b.se.RER=NaN(IRF_length,2);
b_linear.coeff.Y=NaN(IRF_length,1);
b_linear.coeff.G=NaN(IRF_length,1);
b_linear.coeff.RER=NaN(IRF_length,1);
b_linear.se.Y=NaN(IRF_length,1);
b_linear.se.G=NaN(IRF_length,1);
b_linear.se.RER=NaN(IRF_length,1);
var_cell={'G','Y','RER','Y_N'};
var_cell_theoretic={'g_N','y_aggregate','rer','y_N'};

conditioning=true;
% conditioning=false;
if conditioning
    indicator=lagmatrix(cut.U,1)<0.01;
    ending_iter=1;
else
    indicator=true(size(cut.U));
    ending_iter=1;
end

shock_size(1)=1/(loaded_IRFs.(ending_cell{1,ending_iter}).IRFs(loaded_IRFs.(ending_cell{1,ending_iter}).pos.g_N,1,1)./loaded_IRFs.peg_ergodic.Y_base_mean(loaded_IRFs.peg_ergodic.pos.g_N));
shock_size(2)=abs(1/(loaded_IRFs.(ending_cell{1,ending_iter}).IRFs(loaded_IRFs.(ending_cell{1,ending_iter}).pos.g_N,1,2)./loaded_IRFs.peg_ergodic.Y_base_mean(loaded_IRFs.peg_ergodic.pos.g_N)));
scaling_factor(:,1)=shock_size(1)./[loaded_IRFs.peg_ergodic.Y_base_mean(loaded_IRFs.peg_ergodic.pos.g_N)
    loaded_IRFs.(ending_cell{1,ending_iter}).Y_base_mean(loaded_IRFs.(ending_cell{1,ending_iter}).pos.y_aggregate)
    loaded_IRFs.(ending_cell{1,ending_iter}).Y_base_mean(loaded_IRFs.(ending_cell{1,ending_iter}).pos.rer)
    loaded_IRFs.(ending_cell{1,ending_iter}).Y_base_mean(loaded_IRFs.(ending_cell{1,ending_iter}).pos.y_N)];
scaling_factor(:,2)=shock_size(2)./[loaded_IRFs.peg_ergodic.Y_base_mean(loaded_IRFs.peg_ergodic.pos.g_N)
    loaded_IRFs.(ending_cell{1,ending_iter}).Y_base_mean(loaded_IRFs.(ending_cell{1,ending_iter}).pos.y_aggregate)
    loaded_IRFs.(ending_cell{1,ending_iter}).Y_base_mean(loaded_IRFs.(ending_cell{1,ending_iter}).pos.rer)
    loaded_IRFs.(ending_cell{1,ending_iter}).Y_base_mean(loaded_IRFs.(ending_cell{1,ending_iter}).pos.y_N)];

maxNWLag = 5;%floor(4*(length(indicator)/100)^(2/9)); % for Newey-West

for iter=1:IRF_length
    for var_iter=1:length(var_cell)
        y=(cut.(var_cell{var_iter})(iter:end));
        mean_y=nanmean(cut.(var_cell{var_iter})(iter:end));
        X=[ones(length(cut.(var_cell{var_iter})(iter:end)),1), G_positive(1:end-(iter-1)), G_negative(1:end-(iter-1))];       
        [~, se, coeff] = hac(X(indicator(1:end-(iter-1)),2:end), y(indicator(1:end-(iter-1))),'bandwidth', maxNWLag+1,'display','off'); % Newey-West regression
        b.coeff.(var_cell{var_iter})(iter,:)=coeff(2:3)*mean(G)/mean_y;% transform to percent deviations
        b.se.(var_cell{var_iter})(iter,:)=se(2:3)*mean(G)/mean_y;% transform to percent deviations
        X=[ones(length(cut.(var_cell{var_iter})(iter:end)),1), G_shock(1:end-(iter-1))];
        [~, se, coeff] = hac(X(indicator(1:end-(iter-1)),2:end), y(indicator(1:end-(iter-1))),'bandwidth', maxNWLag+1,'display','off'); % Newey-West regression
        b_linear.coeff.(var_cell{var_iter})(iter,:)=coeff(2)*mean(G)/mean_y;% transform to percent deviations
        b_linear.se.(var_cell{var_iter})(iter,:)=se(2)*mean(G)/mean_y;% transform to percent deviations
    end
end
for var_iter=1:length(var_cell)
    b.coeff.(var_cell{var_iter})=b.coeff.(var_cell{var_iter})./repmat(b.coeff.G(1,:),IRF_length,1);
    b.se.(var_cell{var_iter})=b.se.(var_cell{var_iter})./repmat(b.coeff.G(1,:),IRF_length,1);
    b_linear.coeff.(var_cell{var_iter})=b_linear.coeff.(var_cell{var_iter})./repmat(b_linear.coeff.G(1,:),IRF_length,1);
    b_linear.se.(var_cell{var_iter})=b_linear.se.(var_cell{var_iter})./repmat(b_linear.coeff.G(1,:),IRF_length,1);
end

confidence_90  = norminv(0.9 + (1 - 0.9) / 2, 0, 1);
confidence_68  = norminv(0.68 + (1 - 0.68) / 2, 0, 1);

fontsize = 6;
dependent_var_plottitles=var_cell;
dependent_var_ylabels={'percent','percent';'percent','percent';'percent','percent';'percent','percent'};

main_fig=figure;
for shock_iter = 1:2
    var_iter=0;
    for var_pos= [1,3,4]
        var_iter=var_iter+1;
        if shock_iter == 1
            scaling = -1;
        elseif shock_iter == 2
            scaling = 1;
        else
            error('not defined')
        end
        subplot(3,3,3*(var_iter-1)+shock_iter)

        theta_mat=scaling*b.coeff.(var_cell{var_pos})(:,3-shock_iter);
        se_mat=b.se.(var_cell{var_pos})(:,3-shock_iter);
        ci1a_90=0;
        ci2a_90=0;
        ci1b_90 = (theta_mat+confidence_90*se_mat);
        ci2b_90 = (theta_mat-confidence_90*se_mat);
        ci1a_68=0;
        ci2a_68=0;
        ci1b_68 = (theta_mat+confidence_68*se_mat);
        ci2b_68 = (theta_mat-confidence_68*se_mat);

        topa_90 = max(ci1a_90,ci2a_90);
        bottoma_90 = min(ci1a_90,ci2a_90);
        topb_90 = max(ci1b_90,ci2b_90);
        bottomb_90 = min(ci1b_90,ci2b_90);
        topa_68 = max(ci1a_68,ci2a_68);
        bottoma_68 = min(ci1a_68,ci2a_68);
        topb_68 = max(ci1b_68,ci2b_68);
        bottomb_68 = min(ci1b_68,ci2b_68);

        ha1_90 = area(0:IRF_length-1,[bottomb_90 topb_90-bottomb_90],'FaceColor',[204/255 229/255 1],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_90, 'LineStyle', '-')
        hold on
        ha1_68 = area(0:IRF_length-1,[bottomb_68, topb_68-bottomb_68],'FaceColor',[153/255 204/255 1],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_68(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_68, 'LineStyle', '-')
        hold on

        fs=plot(0:IRF_length-1,theta_mat,'b-', 'LineWidth', 2.5);

        theta_mat=scaling*b_linear.coeff.(var_cell{var_pos})(:);
        se_mat=b_linear.se.(var_cell{var_pos})(:);
        ci1a_90=0;
        ci2a_90=0;
        ci1b_90 = (theta_mat+confidence_90*se_mat);
        ci2b_90 = (theta_mat-confidence_90*se_mat);
        ci1a_68=0;
        ci2a_68=0;
        ci1b_68 = (theta_mat+confidence_68*se_mat);
        ci2b_68 = (theta_mat-confidence_68*se_mat);

        topa_90 = max(ci1a_90,ci2a_90);
        bottoma_90 = min(ci1a_90,ci2a_90);
        topb_90 = max(ci1b_90,ci2b_90);
        bottomb_90 = min(ci1b_90,ci2b_90);
        topa_68 = max(ci1a_68,ci2a_68);
        bottoma_68 = min(ci1a_68,ci2a_68);
        topb_68 = max(ci1b_68,ci2b_68);
        bottomb_68 = min(ci1b_68,ci2b_68);

        ha2_90 = area(0:IRF_length-1,[bottomb_90 topb_90-bottomb_90],'FaceColor',[1 0.8 0.8],'EdgeColor','none','ShowBaseLine','off');
        set(ha2_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha2_90, 'LineStyle', '-')
        ha2_90(2).FaceAlpha = 0.5;
        hold on
        ha2_68 = area(0:IRF_length-1,[bottomb_68, topb_68-bottomb_68],'FaceColor',[1 0.6 0.6],'EdgeColor','none','ShowBaseLine','off');
        set(ha2_68(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha2_68, 'LineStyle', '-')
        ha2_68(2).FaceAlpha = 0.5;
        hold on

        fs=plot(0:IRF_length-1,theta_mat,'r--', 'LineWidth', 2.5);

        fs=plot(0:IRF_length-1,scaling_factor(var_pos,3-shock_iter)*loaded_IRFs.(ending_cell{1,ending_iter}).IRFs(loaded_IRFs.(ending_cell{1,ending_iter}).pos.(var_cell_theoretic{var_pos}),1:IRF_length,3-shock_iter),'k-o', 'LineWidth', 2.5);
        fs.MarkerSize=3;
        
        hline(0,'k:')
        if var_iter==2
            if shock_iter == 1
                ylim([0.05 0.35])
            elseif shock_iter == 2
                ylim([-0.5 -0.2])
            end
        end
        xlim([0 IRF_length-1])
        box on;set(gca,'xTick',0:IRF_length-1,'Layer','top','FontSize',fontsize);

        title(dependent_var_plottitles(var_pos),'FontSize',fontsize)
        ylabel(dependent_var_ylabels(var_pos),'FontSize',fontsize)
        xlabel('years','FontSize',fontsize)
    end
end
set(findall(main_fig,'-property','ShowBaseLine'),'ShowBaseLine','off')
print('Figures/MC_IRF_peg_CI','-dpdf')
saveas(main_fig,'Figures/MC_IRF_peg_CI');
save_string='Figures/MC_IRF_peg_CI';

name_vec = {'hike_Y_N', 'hike_FX', 'hike_G', 'cut_Y_N','cut_FX', 'cut_G'};
       
axesH = findobj(main_fig,'Type','axes','-not','Tag','legend'); % This will change the order from which the axes will be saved. If you need to save in the correct order, you will need access to the subplot handles
nAxes = numel(axesH);
newFig = figure;
for k=1:nAxes
    newAxes=copyobj(axesH(k),newFig);
    % Expand subplot to occupy the hole figure:
    [newFig,newAxes]=crop_image(newAxes,newFig);
    print(newFig, [save_string ,name_vec{k}, '.pdf'],'-dpdf')
    delete(newAxes);
end
delete(newFig);


function [newFig,newAxes]=crop_image(newAxes,newFig)
    set(newAxes,'OuterPosition',[0 0 1 1]);
    set(newAxes,'FontSize', 20);
    set(newAxes,'Ylabel',[]);
    set(newAxes,'Xlabel',[]);
    set(newAxes,'Title',[]);
    set(newAxes,'LineWidth', 2);
    set(newAxes,'TickLabelInterpreter','latex')
    tightInset=get(newAxes,'TightInset');
    set(newAxes,'Position',[tightInset(1:2) [1 1]-(tightInset(3:4)+tightInset(1:2))])
    set(findall(newFig,'-property','ShowBaseLine'),'ShowBaseLine','off')
end